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Abstract 

We apply Lindstedt's method and multiple scale per- 
turbation theory to analyze spatio-temporal struc- 
tures in nonlinear Schrodinger equations and thereby 
study the dynamics of quasi-one-dimensional Bose- 
Einstein condensates with mean- field interactions. 
We determine the dependence of the amplitude of 
modulated amplitude waves on their wave number. 
We also explore the band structure of Bose-Einstein 
condensates in detail using Hamiltonian perturbation 
theory and supporting numerical simulations. 

PACS: 05.45.-a, 03.75.Lm,05.30.Jp, 05.45.Ac 
Keywords: nonlinear dynamics, Bose-Einstein conden- 
sates, chaos 

Bose-Einstein condensates (BECs) were ob- 
served experimentally in 1995 using dilute va- 
pors of sodium and rubidium. The macro- 
scopic behavior of BECs at zero temperature 
is modeled by the nonlinear Schrodinger equa- 
tion in the presence of an external potential. 
This model has proven to be an excellent one 
for most experiments on BECs. When the ex- 
ternal potential is spatially periodic (e.g., due 
to an optical lattice, which may be created 
using counter-propagating laser beams), the 
spectrum of the BEC exhibits a band stuc- 
ture (spatial resonance structure). This paper 
utilizes Hamiltonian perturbation theory and 



supporting numerical simulations to study this 
structure in detail. 

1 Introduction 

At low temperatures, particles in a dilute gas can 
reside in the same quantum (ground) state, form- 
ing a Bose-Einstein condensate. 15 ' 23 ' 33 ' 43 This was 
first observed experimentally in 1995 with vapors 
of rubidium and sodium. 4 ' 24 In these experiments, 
atoms were confined in magnetic traps, evaporatively 
cooled to tempuratures on the order of fractions of 
microkelvins, left to expand by switching off the con- 
fining trap, and subsequently imaged with optical 
methods. 23 A sharp peak in the velocity distribution 
was observed below a critical temperature, incidating 
that Bose-Einstein condensation had occurred. 

BECs are inhomogeneous, so condensation can be 
observed in both momentum and coordinate space. 
The number of condensed atoms N ranges from sev- 
eral thousand to several million. Confining traps 
are usually approximated well by harmonic poten- 
tials. There are two characteristic length scales: 
the harmonic oscillator length a^o — y/h/ (mwho) 
[which is on the order of a few microns], where 
uJho = (wiWjWz) 1 ' 3 is the geometric mean of the trap- 
ping frequencies, and the mean healing length x — 
, where n is the mean density and a, the 
(two-body) s-wave scattering length, is determined 
by the atomic species of the condensate. 6, 23,34,43 In- 



teractions between atoms are repulsive when a > 
and attractive when a < 0. For a dilute ideal gas, 
a 0. The length scales in BECs should be con- 
trasted with those in systems like superfluid helium, 
in which the effects of inhomogeneity occur on a mi- 
croscopic scale fixed by the interatomic distance. 23 

If considering only two-body, mean-field interac- 
tions, a dilute Bose-Einstcin gas can be modeled us- 
ing a cubic nonlinear Schrodinger equation (NLS) 
with an external potential, which is also known as 
the Gross-Pitaevskii (GP) equation. BECs are mod- 
eled in the quasi-one-dimensional (quasi-lD) regime 
when the transverse dimensions of the condensate are 
on the order of its healing length and its longitu- 
dinal dimension is much larger than its transverse 
ones. 9-11:23 In the quasi-lD regime, one employs 
the ID limit of a 3D mean-field theory rather than a 
true ID mean-field theory, which would be appropri- 
ate were the tranverse dimension on the order of the 
atomic interaction length or the atomic size. 7 ' 9-11,48 

When examining only two-body interactions, the 
condensate wavefunction ("order parameter") ip(x, t) 
satisfies a cubic NLS, 

ih^t = -[fi 2 /(2m)]fc + g\i>\ 2 i, + V(x)<P , (1) 

where \ip\ 2 is the number density, V(x) is an ex- 
ternal potential, g = [Airh 2 a / m][l + 0(( 2 )], and 
C = VrV^FkF i s the dilute gas parameter. 6, 23,34 Be- 
cause the scattering length a can be adjusted using 
a magnetic field in the vicinity of a Feshbach reso- 
nance, 27 the contribution of the nonlinearity in ^ is 
tunable. 

Potentials V{x) of interest in the context of BECs 
include harmonic traps, periodic potentials ("stand- 
ing light waves"), and periodically perturbed har- 
monic traps. The existence of quasi- ID cylindrical 
("cigar-shaped") BECs motivate the study of peri- 
odic potentials without a confining trap along the 
dimension of the periodic lattice. 37 Experimentalists 
use a weak harmonic trap on top of the periodic lat- 
tice to prevent the particles from spilling out. To 
achieve condensation, the periodic lattice is typically 
turned on after the trap. If one wishes to include the 
trap in theoretical analyses, V(x) is modeled by 

V{x)^V sm{K{x-x Q )) + V 1 x 2 , (2) 

where k is the lattice wave number, Vq is the height 
of the periodic lattice, and x is the offset of the pe- 
riodic potential. (Note that these three quantities 
can all be tuned experimentally.) The periodic term 
dominates for small x, but the harmonic trap other- 
wise becomes quickly dominant. When V\ <C Vq, the 



potential is dominated by its periodic contribution 
for many (20 or more) periods. 14,18 ' 25 (For example, 
when Vq/Vi = 500, k = 10, and x = 0, the har- 
monic component of V(x) essentially does not con- 
tribute for 10 periods.) In this work, we usually let 
V\ = and focus on periodic potentials. Spatially pe- 
riodic potentials have been employed in experimental 
studies of BECs 3 ' 32 and have also been studied the- 
oretically. 2 ' 8 - n ' 18 ' 20 ' 25 ' 40 ' 49 

When the optical lattice has deep wells (large |Vo|), 
an analytical description of BECs in terms of Wannier 
wave functions can be obtained in the tight-binding 
approximation. 41 The Bose-Hubbard Hamiltonian, 
which is a better description than Q in the tight- 
binding approximation, is derived by expanding the 
field operator in a Wannier basis of localized wave 
functions at each lattice site. This Hamiltonian has 
has three contributions: a kinetic energy term yield- 
ing contributions from tunnelling between adjacent 
wells, an energy offset in each lattice site (due, for 
example, to external confinement), and a potential 
energy term characterized by atom-atom interactions 
(that indicates how much energy it takes to put a sec- 
ond atom into a lattice site that already has one atom 
present). One can use the Bose-Hubbard Hamilto- 
nian to examine transitions between superfluidity and 
Mott insulation. 30 

In the present paper, we examine in detail the band 
structure of BECs in shallow periodic lattices us- 
ing Hamiltonian perturbation theory and supporting 
numerical simulations. 44 Our methodology, which 
yields analytical expressions describing the features 
of BEC resonance bands, exploits the elliptic func- 
tion solutions of the NLS in the absence of a poten- 
tial. Note, however, that this paper does not explore 
the chaotic dynamics of BECs. 

2 Coherent Structures 

We consider uniformly propagating coherent struc- 
tures with the ansatz ip(x — vt, t) = R(x — 
vt) exp (i [9(x — vt) — /it]), where R = is the mag- 
nitude (amplitude) of the wave function, v is the ve- 
locity of the coherent structure, 9(x) determines its 
phase, v*o <x V# is the particle velocity, and /i is the 
chemical potential (which can be termed an angu- 
lar frequency from a dynamical systems perspective). 
Considering a coordinate system that travels with 
speed v (by defining x' — x — vt and relabeling x 1 
as x) yields 

i)(x, t) = R(x) exp (i [d(x) -fit]) . (3) 



2 



[From a physical perspective, we consider the case 
v = 0, as V(x') = V(x - vt).] When the (tem- 
porally periodic) coherent structure © is also spa- 
tially periodic, it is called a modulated amplitude wave 
(MAW). 12 ' 13 The orbital stability of MAWS for the 
cubic NLS with elliptic potentials has been studied 
by Bronski and co-authors. 9-11 To obtain stability 
information about the sinusoidal potentials we con- 
sider, one takes the limit as the elliptic modulus k 
approaches zero. 35 ' 45 

When V(x) is periodic, the resulting MAWs gener- 
alize the Bloch modes that occur in the theory of lin- 
ear systems with periodic potentials, as one is consid- 
ering a nonlinear Floquet-Bloch theory rather than a 
linear one. 5 ' 8 ' 20,37 ' 46 In this paper, we employ phase 
space methods and perturbation theory to examine 
MAWs and their concomitant band structure. 

The novelty of our work lies in its illumination 
of BEC band structure through the use of pertur- 
bation theory and supporting numerical simulations 
to examine 2m' : 1 spatial subharmonic resonances 
in BECs in periodic lattices. Such resonances corre- 
spond to spatially periodic solutions ip of period 2m' 
and generalize the 'period doubled' states (in \ip\ 2 ) 
studied by Machholm, et al. 38 which pertain to the 
experiments of Cataliotti, et al.. 19 

Previous theoretical work in this area has focused 
on different aspects of BEC band structure, such 
as loop structure 26,39 ' 52 and hysteresis. 42 In con- 
trast to the coherent structures we consider, these 
authors studied band structure using a Bloch wave 
ansatz. In our notation, they assumed a priori that 
R(x) — R(x + 2tt/h) has the same periodicity of 
the underlying spatial lattice V(x), whereas we have 
made no such assumption and instead use Hamil- 
tonian perturbation theory to study the dynamical 
behavior of R(x). Additionally, the analytical com- 
ponents of these works are confined to two-to-three 

Fourier mode truncations of the Bloch wave dynam- 
^g 26,39,52 

Inserting J2J) into the NLS and equating real 
and imaginary parts yields 



HuRix) = R"(x) 

2m 

tt[0'{x)] 2 +9 R\x) + V{x) 



+ 



R(x), 
(4) 



= \2e'(x)R'(x) + e"(x)R(x)} 

2m 



nonlinear ordinary differential equations 

R' = S, 

" 2 2m/iR 2mg „, 2m 



S' 



,,a R 3 + ^V(x)R. (5) 



i? 3 h 

The parameter c is defined via the relation 

ff (*) = jp> ( 6 ) 

and therefore plays the role of "angular momentum," 
as discussed by Bronski and coauthors. 10 [Equation 
© is a statement of conservation of angular momen- 
tum.] Constant phase solutions, which constitute an 
important special case, satisfy c = 0. 

3 BECs without an External 
Potential 

When V(x) = 0, the dynamical system (0 is 
autonomous and hence integrable, as it is two- 
dimensional. Its equilibria (i?*,^) satisfy S"* = 
and either i?„ = 0, c — or 



gR 6 - h{iR A 



2m 



0. 



(7) 



which can be solved exactly because it is cubic in R 2 . 
When c = 0, one obtains i?* = zty^h/i/g. One thus 
obtains equilibria at (i?*, 0) ^ (0, 0) for g > if fi > 
and g < if /i < 0. 

The eigenvalues of the equilibrium (i?„, 0) satisfy 



A 2 



3c 2 2m/i 6mg 2 



(8) 



When c = and i?* = 0, one obtains A = 
iy— 2rtT4jjH. Additionally, one obtains a center at 
(0, 0) when fi > and a saddle when fi < 0. One also 
obtains saddles at the i?* ^ equilibria for g > 
when /i > and centers at those same locations for 
g < when fi < 0. These latter equilibria are sur- 
rounded by periodic orbits that satisfy R ^ 0. The 
possible qualitative dynamics (for c = 0) are illus- 
trated in Figure n and summarized in Tabled 

To study the dependence of the wave number of 
periodic orbits (centered at the origin) of 101 on the 
amplitude R when V(x) = 0, we employ Lindstedt's 
method and consider null angular-momentum wave 
functions for the case /i > 0. We also assume g = eg, 
where e <C 1 and g = 0(1). The period of R(x) is 
given by 



which gives the following two-dimensional system of 



T 



27r „ 
— = 2tt 

a 



1 



3ffA 2 



(9) 
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where = Rq(£) + 0(e), £ := ax, a = 1 + eai + 
C(e 2 ) is the wave number, A := Rq(0), and 



i? (£) = ,4 cos 



' 2m fi 



(10) 



Note that all periodic orbits are centered at the 
origin when n > 0. When g < 0, the spatial period 
becomes smaller with increasing A. When g > 0, the 
period becomes larger with increasing A. In the latter 
case, the wave number-amplitude relation holds only 
for solutions inside the separatrix, as the trajectories 
are unbounded outside the separatrix and hence not 
periodic. 

Before deriving the wave number-amplitude rela- 
tions when V(x) ^ 0, we comment briefly on the pro- 
ceeding results. The spatial period for small c ^ 
is similar to (JjJJ, but it cannot be estimated as easily 
because equation J5J now includes a term of order 
0(R~ 3 ) with coefficient c. Although © can be com- 
puted exactly in terms of elliptic functions, here we 
are interested in elucidating the qualitative dynamics 
of the MAWs of interest as well as establishing the 
methodology to be employed in the presence of poten- 
tials V(x). We will utilize elliptic function solutions 
in Section|S|in our detailed study of band structure. 44 
The physical relevance of elliptic functions to BECs 
has been discussed by Carr and collaborators. 16 ' 17 

4 BECs in a Periodic Lattice 

To study the wave number-amplitude relations of pe- 
riodic orbits in the presence of external potentials, 
we expand the spatial variable x in multiple scales. 
We define "stretched space" £ := ax as in the inte- 
grable situation and "slow space" rj := ex. We con- 
sider potentials of the form V(x) — eV(^,rf), where 
V(£, rj) = V sin [k(£ - &)] + V\{rj) and Vx, which is 
of order 0(1), is arbitrary but slowly varying. Cases 
of particular interest include Vx = (periodic poten- 
tial) and Vx = Vx(rj — t]q) 2 (superposition of periodic 
and harmonic potentials). 

When k 7^ ±2(3 :— ±2^j2mfi/h, the equations of 
motion for the slow dynamics of j5J| with c = are 

" A '-^ B -B Bc2 -m m " } - 



AVUl), (11) 



di] 

dB — _ p a j± _|_ ^Pg ^2 _|_ 
dr/ 8/ih 2fih 

where g = eg. The leading-order expression for the 
amplitude is 

Ro(Z, rj) = A( V ) cos(/3£) + B( v ) sin(/3£) , (12) 



where C 2 :— A(r]) 2 + B(rj) 2 is a constant. The dy- 
namical system (I11H is autonomous when V\(x) = 0. 
Equilibria (A*,B*) ^ (0, 0) of ((TT|) correspond to pe- 
riodic orbits of JSJ) with c = 0. The equilibrium value 
of the squared amplitude is denoted C 2 — A 2 + B 2 . 

Converting to polar coordinates with A(rf) = 
C cos(<j)(r])) and B(rj) = Csm.(<p(ij)) and integrating 
the resulting equation yields 



4>(v) = 0(0) 



V 



(13) 

The wave number of the periodic motion is given 



by 

a ^ = 1 -ik c2 -^h v ^ +0 ^- (14) 

When k = ±2/3, we show that the slow flow equa- 
tions have an extra term due to resonance. With- 
out loss of generality, we let k = +2/3, as changing 
the sign of Vq produces the k = —2/3 case. When 
m = 0.5, /i = 10, and h = 1, for example, one obtains 
this resonant situation for k = ±10. Additionally, we 
show that C is no longer constant in this resonant 
situation. 

In polar coordinates, the slow flow equations are 



=ai(j+ 8Hh C + 4^ Sm[2( ^ 



/%>)] 



C = --^Ccos[2(0- 



mo 



(15) 



Integrating the equation for C" yields 



C = Cq exp 



VoP 



cos[2(0(r ? ) - Ko)]dr] 



(16) 



which one may then insert into the equation for the 
angular dynamics. 

To determine equilibria, one puts C = <f>' = 0. 
From C = 0, one determines that equilibria (C*, 0*) 
satisfy 



(2j + 1)tt 



■/%, je {0,1,2,3}, 



(17) 



which is independent of the scattering coefficient. In- 
serting lfT7|) into C — yields the wave number- 
amplitude relation 



a R (C) 



?s~t\ Vq 



0(e 2 ) 



(18) 



for periodic orbits of JSJ. In (|18fl. the minus sign 
is obtained when j £ {0, 2}, and the plus sign is ob- 
tained when j G {1,3}. Equation QlSj l is valid for 2 : 1 
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spatial resonances. We examine 2m' : 1 resonances for 
integer m' in Section [5] using Hamiltonian perturba- 
tion theory and the elliptic function solutions of JSJ 
when Vq = 0. 

To examine the spatial stability (i.e., stability with 
respect to spatial evolution) of these periodic orbits in 
the presence of resonant periodic potentials, we com- 
pute the spatial stability of equilibria of l|15f) when 
V\(x) = 0. The eigenvalues of the periodic orbits are 



(19) 



We show numerical simulations for (JjJ in the presence 
of a periodic potential in Figure In this situation, 
(0 is a nonlinear Mathieu equation. 45 ' 46 ' 53 FigureEJl 
shows the coherent structure for the trajectory with 
K = 100 and Vq — 10. Figure |3] depicts a Poincare 
section describing the dynamics of S5 Rb, for which 
a = —0.9 nm. Figure ^ depicts spatial profiles of 
the coherent structures corresponding to the locally 
chaotic and globally chaotic trajectories in Figure 

5 Subharmonic Resonances 

In this section, we analyze spatial subharmonic res- 
onances and the band structure of repulsive BECs 
with a positive chemical potential. We perturb off 
the elliptic function solutions of the underlying inte- 
grable system in order to study 2m' : 1 spatial res- 
onances with a leading-order perturbation method. 
Perturbing off simple harmonic functions, by con- 
trast, requires a perturbative method of order m' 
to study 2m' : 1 resonances. At the center of the 
KAM islands, we observe 'period-multiplied' states. 
When m! = 1, one obtains period-doubled states in 
ip. As verified numerically in Section [H] our qualita- 
tive results are excellent. Given that our method is 
a leading-order one, our quantitative results are also 
remarkably good. 

Recent work by Machholm and coauthors 38 on 
period-doubled states (in \ip\ 2 ) follows up experi- 
mental studies by Cataliotii and coauthors, 19 who 
observed superfluid current disruption in chains of 
weakly coupled BECs, which is related to the de- 
pendence of the dynamical instability of Bloch states 
on the magnitude of particle interactions. Period- 
doubled states, which may be interpreted as soliton 
trains, arise from dynamical instabilities of the en- 
ergy bands associated with Bloch states. 38 In the 
present work, we offer a dynamical systems perspec- 
tive on period-doubled states and their generaliza- 



tions. Our theoretical and computational analysis re- 
veals period- multiplied solutions of the GP QJ. The 
existence of these wave functions can be explored ex- 
perimentally. 

A detailed examination of the band structure of 
BECs in periodic lattices requires a more intricate 
perturbative analysis than that discussed earlier in 
this work. Previous authors have concentrated on 
numerical studies of band structure. 8 ' 20 ' 37 The ap- 
proach we take, on the other hand, is to analyze 
the spatial resonance structure that arises from the 
nonlinear Mathieu equation obtained upon the ap- 
plication of a coherent structure ansatz to the cubic 
NLS. We examine situations with null angular mo- 
mentum (c = 0), but one observes similar behavior 
when c ^ when R is away from the origin. The ana- 
lytical approach we employ was introduced by Zounes 
and Rand 53 for g < and /i > (see Figure H 3 ) , the 
technically easiest case to consider. Their study of 
nonlinear Mathieu equations is directly applicable to 
BECs. Our work is an extension of their work to the 
situation g > 0, fi > (see Figure^), the second 
easiest case to consider. We study this case in detail 
and also apply the results of Zounes and Rand to at- 
tractive BECs with a positive chemical potential. We 
briefly discuss attractive BECs with a negative chem- 
ical potential (see Figure^), the technically hardest 
case to consider. Note that this paper does not ex- 
plore the chaotic dynamics of BECs, which is an im- 
portant open issue. 31, 45,46,51 

Let xq = — 7t/(2k) and V\{x) = so that 



V{x) = Vq COs(kx) . 



(20) 



When c = 0, the equations of motion J2J), JSJ for the 
amplitude of the coherent structure © take the form 



R" + SR + aR 3 + eR cos(kx) 



where 



2mfx 



2mg 



2m, 
If 



■Vo- 



(21) 



(22) 



(Note that the perturbation parameter e is not the 
same as the parameter e employed earlier.) The pa- 
rameters n, Vo, K, and a (and hence g) can all be 
adjusted experimentally. When e = 0, solutions of 
(|21|) can be written exactly in terms of elliptic func- 

t j ons .21,22,35,45,50,53 



R = apcn(u, k) . 



(23) 
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where 



u = U\X + Uq , Uj = 8 + ap , 

2(5 + ap 2 ) ' 
ux > 0, p>0, fc 2 eM, ere {-1,1}, 



(24) 



and Mo is obtained from an initial condition (and can 
be set to without loss of generality). We consider 
u\ 6 1 in order to study periodic solutions. One 
can use argument transformations to study solutions 
with complex u\. When k 2 £ (1, oo), one makes sense 
of the cn function with a reciprocal modulus trans- 
formation. 21 ' 22 When k 2 < 0, one employs a recip- 
rocal complementary modulus transformation, which 
we discuss below. 

Equation (|21|1 can be integrated when e — to 
yield the Hamiltonian 

ii?' 2 + ^SR 2 + jaB 4 = h , (25) 

with given energy h. With l|24|) . one computes 

^ 2 ^ + ^ = ?(I^< ^ 

where k' 2 :— 1 — k 2 . Earlier in this paper, we enu- 
merated the different possibilities for the qualitative 
dynamics of (|21|l in terms of the signs of p and g (and 
hence in terms of the signs of S and a) . 

5.1 Repulsive BECs with a Positive 
Chemical Potential 

We first consider in detail the case g > 0, fi > 0, for 
which S > 0, a < 0. For notational convenience, we 
sometimes utilize a' :— —a. This analysis involves 
a considerable amount of elliptic-function manipula- 
tion, but we are rewarded in the end by a much more 
effective perturbation theory than can be obtained 
by employing trigonometric functions. 

The center a t (0, 0) satisfies h = p 2 = k 2 = 0. The 
saddles at (zty^S/a', 0) and their adjoining separatrix 
satisfy 



4a 



R 



(27) 



The sign a = +1 is used for the right saddle, and a = 
— 1 is used for the left one. Within the separatrix, all 
orbits are periodic and the value of a is immaterial. 



5.1.1 Action-Angle Variable Description and 
Transformations 

For this choice of parameters, k 2 £ [— oo, 0], so elliptic 
functions are defined through the reciprocal comple- 
mentary modulus transformation, 21 ' 22 which relates 
the (it, k) coordinate system to another coordinate 
system, which we denote (w, fo)- To tranform be- 
tween these two coordinate systems, one uses the fol- 
lowing relations: 



cn(tt, k) — cd(u>, kz) , 
dn(tt, k) = nd(w, fo) , 
sn(u, k) = fc 2 sd(ui, £2) , 

= — ii = hint K = hL K n K = 1 

h 



k' 



u = k' 2 w, K= k' 2 K 2 , E = -tjEi . 

2 ^2 



(28) 



Here, K = K(k) denotes the complete elliptic inte- 
gral of the first kind, E = E(k) denotes the complete 
elliptic integral of the second kind, and items with 
the subscript '2' denote the analogous quantities in 
the (w,k2) coordinate system. 21 ' 35 ' 50 

We rescale l|21|l using the coordinate transforma- 
tion 

j r 

X = VSx, r =\—R (29) 



to obtain 







(30) 



when V(x) = 0. (Note that in this analysis, the quan- 
tity x does not represent the mean healing length.) 
In terms of the original coordinates, 



Imp 



(31) 



The rescaling applied for other choices of S and a 
differ slightly from that in l|29Jl , so that the arguments 
of their associated square roots are positive. 
The Hamiltonian corresponding to l|3UI) is 

H (r, s) = X -s 2 + l -v 2 - ir 4 - h , h G [0, 1/4] , (32) 

where s := r' = dr/dx- Additionally, p 2 6 [0,1], 
k 2 S [0, 1] (corresponding to k 2 e (— oo, 0] in the 
original coordinates), and 



(33) 



2(p 2 1) ' 

With the initial condition r(0) = p, s(0) = 0, which 
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implies that uq — 0, solutions to (|30[) are given by with the corresponding Hamiltonian 

r( X ) = pen ([1 - p 2 ] 1/2 X , *) , S ' X) = ^ o(r ' S) + £i?l(r ' S ' X) 

. 2 1 1/2 «n r n - « 2 1 1/2 v fc^ = I s " + \r 2 ~r 4 + ±t* cos ( ^= X 



The period of a given periodic orbit T is 



xdn([l-p 2 ] 1/2 X ,fc). (34) (42) 

In action-angle coordinates, this becomes 

H{^J,x) = h(J)+€h 1 (9,J, X ) 
T(k) =£d X = , (35) = l - P {jf Ip( J) 4 (43) 

where AK{k) is the period in u of cn(«, fc). 50 The + 2S P ^ ^ WW 71 "' k ) cos 

frequency of this orbit is 

One obtains a second action-angle pair (0, j) us- 



7ryT — p2 i n g the canonical transformation (<£>, J) — ► (</>, j) 

n ( k ) = 2K(k) ' ^ 

defined by the relations 

Let Th denote the periodic orbit with energy h = j(J) = -p(J) 2 , $(0, j) — — — , (44) 



Ho(r,s). The area of phase space enclosed by this 

orbit is constant with respect to X , so one may define wnere 

the action 21 - 22 ' 29 - 45 



2- ' ' — J'(j) 



2tt J r 2tt 



fc 2 = 



J:=±I sdr=±- ( T{k) [s(x)?d X , (37) 2 J 1 



J(j) = ^VT^2j [E(j) - (1 - j)#(j) 



which in this case can be evaluated exactly: ~ . 2 r . „ , 2 „ r , . , , 

K(j) = -K[k(j)}, E(j) = -E[k(j)}. (45) 



7T 7T 



J = 4v/ ^ ^ [E(k) - (1 - p 2 /2) if (fc)] . (38) Additionally, 
The associated angle 29, 31:36,51 in the canonical trans- J'(j) := — = \J\ — 2jK(J) — ^ . (46) 



formation (r, s) ► (J,^) is 

$ := $(0) + 0(fc)X- (39) 



Note that J ~ j for small-amplitude motion. Fur- 
thermore, j = at the origin, and j = 1/2 on the 
„ separatrix. 

The frequency monotonically decreases as fc The Hamiltonian g3J becomes 

goes from — oo to [that is, as one goes from the 

separatrix to the center at (r,s) = (0,0)]. With this H((f>,j,x)=H Q (j) + eH 1 (<f>,j,x) (47) 
transformation, equation (13411 becomes / _ \ 

= 7 - 7 + - ;cn 2 — fc cos — v 
r(J, $) = p(J) cn (2X(fc)$/7r, fc) , <T U'(j) / Vv^ 



«(X) = -p(J)Vl-P^) 2 sn(2^(fc)$/7r,fc) 



Because we have used elliptic functions rather than 



x dn (2if(fc)$/7r, fc) , (40) trigonometric functions, all results are exact thus 

far. 53 

where fc = k(J). 

After rescaling, the equations of motion for the Perturbative Analysis 



forced system (|21|) take the form 



A subsequent 0(e) analysis at this stage allows one to 



. e ( k \ / * ~t \ study 2m' : 1 subharmonic resonances for all ml G Z. 



r r 7 ^ S COS \y/S / r ^ con t ras t, had we undertaken this procedure with 
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trigonometric functions (which would have entailed 
a perturbative approach from the beginning), an 
<D(e m ) analysis would be required to study 2m! : 1 
subharmonic resonances of (|2I[1 . 

The Fourier expansion of cn is given by 



cn(u, k) = 



2tt 



kK{k) 



cos 



(2n + l) 



2K(k) 



where the Fourier coefficients b n (k) are 
1 



(48) 



b n (k) = -sech[{n + l/2)TrK'(k)/K(k)] , (49) 



and K'(k) := K(\/l — k 2 ) denotes the complemen- 
tary complete elliptic integral of the first kind. 1, 50, 53 
In the present situation, 



(2n + l) 



J'(J) 



where 




wm Km] 



= B {j) + J2Bi cos 



(50) 
(51) 



2l<f) 



(52) 

where the Fourier coefficients Bi(j) are obtained by 
convolving the previous Fourier coefficients (|5I|) with 
each other. 53 

Before proceeding, it is important to discuss the 
computation of the coefficients Bi(j), which require 
some care. Using the Elliptic Nome 1 ' 54 



q(k) := e -K'{k)/K(k) : 
the Fourier coefficient (|49|l is expressed as 
b n (k) 



g(fc)n+l/2 + g ( fc )-(n+l/2) 



(53) 



(54) 



One then expands Bi(j) in Taylor series about j = 0. 
In this computation, one finds that the coefficients of 
even powers of j in Bi(j) are the same as when g < 0, 
[i > and that odd powers have the opposite sign. 
This distinction lies at the root of the qualitatively 
different dynamics in the two cases, which we will 
discuss in Section Xb.2\ Recall that their underlying 
integrable dynamics are depicted in Figure ^ 



After the Fourier expansion, the perturbative term 
in the Hamiltonian (|47|l is 



Hi(<j),j,x) = ^B {j) cos [-^X 



oo 



i=i 
cos 



2l<f) 



21(f) K 

J 7 U) + 7s ) 

K 



(55) 



There are infinitely many (subharmonic) resonance 
bands, 31,45, 51 each of which corresponds to a single 
harmonic in the perturbation series H55|l. To isolate 
individual resonances, we apply a canonical, near- 
identity transformation 29,31 ' 45 ' 51 ' 53 to the Hamilto- 
nian H = Hq + eH±. This transformation is given 

by 



P 



dWi 

BP 
dWi 



0(e 2 ), 



Oie 1 



(56) 



where the generating function W\ is 



Wi 



pBq(p) sip ( k 

p 00 

5 B,(p) 

1=1 l^tm' 



1=1 , lj= 
sin 



2iVsn(p) 



( 2l Q _ 
\J'(P) Vs x 

- 2iVsn(p) 



(57) 



To obtain l|57|) . one uses the fact [from (|46[) ] that 
Q(P) = (1-2P)/J'(P). 

The resulting Hamiltonian is 

K(Q, P, x) = K (P) + eK^Q, P, x ) , 
K Q (P) = P - P 2 = H (P) , 



K 1 (Q,P,x)=H 1 (Q,P,x) + {H ,Wi} 



dWi 

dx 



(58) 



where {Ai,^} denotes the Poisson bracket of A\ 
and Ai. For the present choice of Wi, one obtains 
the resonance Hamiltonian 



K(Q,P, X ;m')=P-P 2 



Y8 PB ^ p)co \np)-7i' 



0(e 2 ). 
(59) 
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The choice of the generating function (|57|) eliminates 
all resonances from the Hamiltonian K except the 
2m! : 1 resonance. In focusing on a single resonance 
band in phase space, one restricts P to a neigh- 
borhood of P m >, which denotes the location of the 
m'th resonant torus (associated with periodic orbits 
in 2m! : 1 spatial resonance with the periodic lattice) . 

5.1.3 Resonance Relations 

Resonant frequencies arise when the denominators of 
the terms in W\ vanish, 31 ' 45, 53 which yields the equa- 
tion 

4= = ±2m'n(P m ) (60) 

for the resonance of order 2m! : 1. As ft < 1 is a 
decreasing function of P e [0,1/2), the resonance 
band associated with 2m! : 1 subharmonic spatial res- 
onances is present when 



V6 



< 2m! . 



(61) 



For example, when k = 2.5 and 8 = 1, there are 
resonances of order 4:1, 6:1, 8:1, etc, but there are 
no resonances or order 2:1. When k = 5 and 8 = 1, 
there are resonances of order 6:1, 8:1, 10:1, etc, 
but there are no resonances of order 2 : 1 or 4 : 1. In 
terms of the original parameters, the condition (|61|l 
describing the onset of 2m! : 1 resonance bands takes 
the form 



k < 2m! 



2m/i 



(62) 



If the lattice V(x) has a smaller wave number 
(larger periodicity), then the chemical potential /i 
must be smaller for a given resonance to occur. As n 
is decreased for a fixed y, (i.e., 8) or as fi is increased 
for a given lattice size k, resonance bands of lower 
order emerge from the origin and propagate in phase 
space. Consequently, a sufficiently high order reso- 
nance is always present in (|21|) , but a given number of 
low-order ones may not be. Lower-order resonances 
occupy larger regions of phase space, so l|61[) also in- 
dicates the volume of phase space affected by spatial 
resonances. We will illustrate this in more detail in 
Section with numerical simulations. 



change of coordinates. 45 ' 53 Toward this end, we de- 
fine the generating function 



F(Q,Y, X ;m') = QY- 



2m' VS 



J(Y)X, (63) 



which yields 



P = ^(Q,Y, X )=Y, 
dF 

Z = w (Q,y,x) 



2m'V8 



J'(Y) X = Q- 



2m'VS 



J\P)X- 
(64) 



(Note that in this analysis, £ does not represent 
stretched space, as it did in our multiple scale expan- 
sion.) With this final transformation, the resonance 
Hamiltonian (|59J) becomes 

dF 

K m , (£, Y) = K(Q, P, X ; m') + — (Q, Y, x ) 



Y — Y — 



2m' VS"' 

+ ^(r)cos(|^), (65) 

which is integrable in the {Y, £) coordinate system. 
In [R, 5)-space, level curves of K m i correspond to 
invariant curves of Poincare sections of (|21|l . which 
are defined by strobing the system when the spatial 
variable takes the values x n = 2nir/K. 

We now provide an analytical description of the 
resonance bands under discussion. In particular, we 
compute the locations and type of equilibria and 
width of resonance bands as functions of the param- 
eters 8, e, and k, and hence of /i, Vo, and k. Such 
bands emerge from the action P = Y = Y m i, which 
designates the location of the m'th resonance torus 
in phase space and is determined by the resonance 
relation Iptljl : 



-= = 2m'n(y m >). 



(66) 



5.1.4 Analytical Description of Resonance 
Bands 

To further examine the resonance structure of l|21|) . 
we make (|59fl autonomous via another canonical 



This resonance band is associated with periodic or- 
bits in 2m! : 1 spatial resonance with the periodic 
lattice. The resonance torus is filled with degenerate 
periodic orbits that split 31 ' 51 into 2m! saddles and 
2m! centers when a perturbation is introduced. 
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From H(j5|) . one obtains Hamilton's equations 
Y' 



8K m , eYB m > (Y) . (2m'£ 
■ sm 



J'{Y) 



dK m 
dY 



= 1-2Y- 



2m'VS 



e 

Yd 



(YB m >(Y))' cos 



J'(Y) 

1W) 



T"(Y) 

2m' £, — -jYB m ' (Y) sin 



Equilibria satisfy either Y = or 
2m'£ 



2m'i 
JUYj 



sm 



J'(Y) 



= 0. 



(67) 



(68) 



They also satisfy 
£' = = 1 - 2Y 



— — = J'(Y)±-[YB m ,(Y)]' , 
2mW5 26 



where the sign ± in l|69|l arises from 
2m'£ 



j'on 



0. 



(69) 



(70) 



Using J'(Y) = VI - 2YK(Y), equation ^ is writ- 
ten 

1-2Y- ^-^VT^2YK(Y) ± ± [YB m > (Y)}' = . 

(71) 

When g < and \x > 0, the + case yields a saddle 
and the — case yields a center. In the present situa- 
tion (g > 0, /i > 0), this holds for odd m'. When m! 
is even, — is a saddle and + is a center. 

At equilibria, the action Y takes the value 

Y e = Y m , + eAY + 0(e 2 ) = Y m . ± 0(e) , (72) 

with the signs as in (|71|) . However, note that Y" c > 
Y m i > Y s , just as for g < 0. One inserts into (|7T|l 
and expands the result in a power series. At order 
0(e°) = 0(1), this reproduces the resonance relation 
(I66|) . At order 0(e), one obtains 



n(y m ,)vT^2i^7A''(y m o - 1 



, (73) 



where saddles Y s use the + sign and centers Y c use 
the — sign when m! is even, and the opposite is true 
when m! is odd. When m' is even, Ay > 0, but 



Ay < when m! is odd. Additionally Y c is always 
larger than Y s (for both signs of e). 

Resonance bands occupy a finite region of phase 
space bounded by a pendulum-like separatrix. When 
a perturbation is introduced, trajectories outside the 
separatrix behave almost as they would in the ab- 
sence of a perturbation, so it is important to estimate 
the width of resonance bands, which emerge at action 
values satisfying the resonance relation l|66l) . Because 
of the direction of the inequality in (|62[1 . this is more 
of a condition for non-existence of given resonances. 
[See the discussion following equation (|6~TJl .] For a 
given set of parameters, there will always be reso- 
nances of sufficiently high order (i.e., for a sufficiently 
large m'). However, as we illustrate numerically be- 
low, there are parameter regions in which no 2 : 1 
resonances exist, regions in which no 2 : 1 or 4 : 1 res- 
onances exist, etc. This behavior contrasts markedly 
with that observed when g < 0. 53 In that situation, 
there exist parameter regions in which only 2 : 1 res- 
onances exist, regions in which only 2 : 1 and 4 : 1 
resonances exist, etc. 

We now show that the width of a resonance band 



O 



$m' C-^m' ) | ^| 



(74) 



for perturbations of size e = — 2mVo/H 2 . 

The separatrix of interest passes through the sad- 
dle point Y s , and the maximum extent of the reso- 
nance band occurs at the same phase £ as the asso- 
ciated center, so 



h',„- ( ens | -J^Jr I L,Y ) . 



J'(Ys. 



K m > cos 



2m'i 
J^Y) 



i,y (75) 



when m' is odd and 

2m'i 



Km' cos 



J'{Y S 



-1, Y = Y f 



K m t cos 



2m'£, 

when m! is even. This implies that 



l,y (76) 



(y - y) - (y - y) 2 - ^^(J(Y) J<X.) 

±^(YB m ,(Y)+Y s B m ,(Y s )) = 0, (77) 

where the + sign holds for odd m! and the — sign 
holds for even m'. (Only the + case needs to be 
considered when g < and /i > 0.) 
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Solutions Y"* of (|77ll are perturbations to Y m i of the 
form 

Y, =Y m . +We' 1 + 0{^) (78) 

for an appropriate choice of 7, to be determined by a 
self-consistency argument. (When e < 0, one writes 
(|78|l with (— e) 7 instead. Everything stated here is 
otherwise the same in that situation.) In this analy- 
sis, one uses the fact that Y s = Y m > ± e(AY~) + 0(e 2 ), 
where the + sign is for odd ml and the — sign is for 
even ml . 

To find W and 7, we insert F* and Y s into (|77|l 
and expand the resulting expression in a power series 
about e = 0. During this process, one obtains 

Y 2 = Y r 2 n , + 2Y m ,W^ + W 2 e 2 ^ , 
J(Y*) = J{Y m .) + e 7 WJ\Y m ,) 

+ e 2 ~<W 2 J"{Y m ,) + 0{e i ~<), (79) 

which shows that that the only suitable value of 7 
is 1/2. Equating terms of order 0(1) yields no new 
information. Equating terms of order 0(e 1 ^ 2 ) yields 
the resonance relation H66JI • Equating terms of order 
0(e) shows that 



W = 



±- 



1 



2m 



\r & J"(Y m >) 



1/2 



(80) 



where the + sign occurs for odd ml and the — sign 
occurs for even ml . Therefore, the miminal action of 
the resonance band is 



Y mm =Y m , -V~eW + 0{e). 



and the maximal action is 

Y max =Y m .+^lW + 0{e). 



(81) 



(82) 



Y„ 



The width of the resonance band is Y max 

In Section we compare these analytical results 
with numerical simulations. 



When 6 > 0, a > 0, the phase space of the integrable 
problem contains no separatrix, and the entire space 
is foliated by periodic orbits. (See Table ^) This 
choice of parameters also leads to the simplest ap- 
plication of the perturbation technique described in 
Section |5~T1 In this case, k S [0, 1], so one need not 
apply a modulus transformation in the elliptic func- 
tion solution. One may also set a = 1. 

We refer the reader to Zounes and Rand 53 for de- 
tails. Here, wc highlight a few results that we wish 
to contrast directly. When g < and \i > 0, the 
resonance relation one obtains is 



-j= = 2m'n a (P m ,), 



(83) 



where the frequency fl a (P) has a similar form to that 
of fl described above. In this situation, £l a (P) > 1, 
so subharmonic periodic orbits are present when 



V5 



> 2m' , 



(84) 



which is the reverse inequality as that derived in the 
repulsive case. Hence, there exist regimes in which 
only 2 : 1 resonances are present, only 2 : 1 and 4 : 1 
resonances are present, etc. In terms of BEC parame- 
ters, the condition (|84|l describing the onset of 2ml : 1 
resonance bands takes the form 



k > 2m! 



2m/j 



(85) 



Because the inequalities in (|62|l and l|85() are oppo- 
sitely directed, adjustments to k and /x have the op- 
posite effect in these two cases. 

Additionally, in this case there is no alternating of 
signs in the location of saddles and centers in res- 
onance bands, as there is when g > and /i > 0. 
Because the attractive case with a positive chemical 
potential is simpler than the one we studied, Zounes 
and Rand 53 were able to obtain better predictions de- 
scribing the location of saddles and centers and the 
width of resonance bands from a perturbation anal- 
ysis like that discussed in Section |5~T1 



5.2 Attractive BECs with a Positive 
Chemical Potential 

Zounes and Rand 53 considered l|21|l when S > and 
a > (in other words, n > and g < 0), which is de- 
picted in Figure ^p. They did not consider the appli- 
cation of their analysis to Bose-Einstein condensates, 
so we presently interpret their results in this new light 
and compare it to our analysis of the repulsive case. 



5.3 Attractive BECs with a Negative 
Chemical Potential 

The most difficult case to consider is that of attrac- 
tive BECs with a negative chemical potential. In 
(|2U, a > and 8 < (i.e., /i < 0), so the inte- 
grable dynamics exhibit two homoclinic orbits. (See 
Figure^;.) The perturbative approach used in this 
paper must be applied separately inside and outside 
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the separatrix. Orbits inside the separatrix satisfy 
h < 0, those on the separatrix satisfy h — 0, and 
those outside the separatrix satisfy h > 0. 

Inside the separatrix, k G (l,oo), so one must ap- 
ply the reciprocal modulus transformation to the ar- 
guments of the elliptic functions i(2"3|) , I^H • The sign 
of a determines whether one is considering perturba- 
tions of periodic orbits in the right half or left half 
of the phase plane. To utilize our perturbative anal- 
ysis outside the separatrix, one must expand elliptic 
functions and elliptic integrals in power series about 
infinity, where k = 0. This requires delicate numeri- 
cal computations of Laurent series coefficients. 

In principle, one can overcome the increased tech- 
nical challenges present in this third case (which is 
also of interest) and apply the same analysis as in Sec- 
tion [O] but the lengthy calculations involved would 
entail a separate publication. 

6 Numerical Simulations 

To compare the analytical results in Section [S] with 
numerical simulations, wc utilize (R, S) coordinates 
with m = 1/2 and h = 1. To lowest order in e, the 
change of variables Y — ► P — ► j is a near-identity 
transformation, so Y = j + O(e). Recall from (|32|l 
that 

■ 1 2 

J= 2 P > 

H = h{j)=j-f = \s 2 + \r 2 -\r\ (86) 

where s = dr/dx- For this comparison, we let a' — 
1 and vary k, 5 = fi, and e = —2mVo/h 2 . Recall 
additionally from Ij29(l that 

/ q „ ct' 1 / q . yV . 

r = \ ir R= vt r > s = ~\Hr R = hr R ■ 

y h[i V o n V 2m o 

(87) 

6.1 Methodology 

Before discussing our results, we briefly overview our 
comparison procedure. 

The "exact" locations of saddles and centers and 
sizes of resonance bands were determined using direct 
numerical simulations of Poincare sections of l|21l) . 
The surface of section we employed satisfies x n = 
Itvk j k (n G Z), which consists of integer multiples of 
the periodicity of the sinusoidal forcing in (|21|l . In our 
simulations, the variable kx is periodic, so the surface 
of section is defined by the condition kx = 0. We used 



this framework to find saddles, centers, and resonance 
band sizes (i.e., separatrix widths) empirically. 

To obtain our predictions, we employed the res- 
onance Hamiltonian l|t)5|) . whose level curves corre- 
spond to invariant curves of Poincare sections. As 
each trajectory yields a level set of this Hamiltonian, 
we solved K' m = constant numerically at appropri- 
ate energy values to obtain predictions for the loca- 
tions of saddles and centers and the size of resonance 
bands; these latter quantities are determined from 
the widths of separatrices in . For these compu- 
tations, we expanded elliptic functions and elliptic in- 
tegrals in Taylor series and subsequently transformed 
these results to (R, 5)-space to compare these calcu- 
lations with our empirical ones. We also predicted 
the locations of saddles and centers I|72I73|I and the 
size of resonances bands <|80I81I82|) using the predic- 
tions obtained from further perturbation expansions. 
We again tranformed back to (i?, S) space to com- 
pare this second set of predictions with our empirical 
results. 

6.2 Primary Resonances 

Our comparison between theory and numerics for pri- 
mary resonances is summarized in Tables [21 and 

Consider first K = 1.5 and 5=1. Poincare sections 
and level sets of the resonance Hamiltonian K\ [in 
units of £/J'(Y)] are depicted in Figure[S] The results 
of our comparison between perturbation theory and 
numerical simulations are summarized in Tabled 

We do relatively well in locating saddles and ex- 
tremely well in locating centers. This is especially 
significant in light of the fact that many canonical 
transformations were required to obtain our analyt- 
ical predictions. Although the requisite calculations 
are complicated, we are rewarded by excellent quali- 
tative agreement and good (and sometimes excellent) 
quantitative agreement. For e = 0.05, the numerical 
resolution of the location of the saddles was prob- 
lematic, so a direct comparison is necessarily less ac- 
curate. As a result, a range of values is sometimes 
indicated for the numerically determined location of 
saddles. Such difficulties with direct numerical simu- 
lation emphasize the importance of using qualitative 
analytical methods to study the features of resonance 
bands. 

Our comparisons between perturbation theory and 
numerical simulations for k = 0.75, 5 = 0.2 are sum- 
marized in Table |21 

If desired, one can improve these quantitative pre- 
dictions by including higher-order contributions in 
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the perturbation expansions. 

6.3 Secondary Resonances 

Our comparison between theory and numerics for sec- 
ondary resonances is summarized in Table 

We study 4 : 1 resonances for k = 2.5 and 5=1. 
No 2 : 1 resonances exist for this choice of parameters. 
The resonance Hamiltonian is depicted for e = 0.05 
in Figure [fj] The corresponding Poincare section is 
shown in Figure [7| 

When e = 0.01, we observe numerically that 
centers are located at approximately (i?, S) — 
(±0.691, ±0.332). With R = ±0.691, we predict a 
value of S = ±0.32384. With S = ±0.332, we pre- 
dict a value of R — ±0.68362. These predictions are 
remarkably good, as we have used leading-order per- 
turbation theory to derive analytical predictions for 
4: 1 (secondary) resonances. However, they are not as 
good as those obtained for the location of saddles in 
this case or the location of centers for 2 : 1 (primary) 
resonances. 

When e = 0.05, numerical simulations sug- 
gest that centers are located at about (R, S) — 
(±0.697, ±0.330). Using R = ±0.697 leads to a pre- 
diction of S = ±0.31912. Using S = ±0.330 leads to 
a prediction of R = ±0.68721. 

6.4 Tertiary Resonances 

Our comparison between theory and numerics for ter- 
tiary resonances is summarized in Table 

We consider 6 : 1 resonances for k = 3.8, S = 1, and 
e = 0.01. No 2 : 1 resonances exist for this choice of 
parameters, but 4 : 1 resonances do exist. The reso- 
nance Hamiltonian is depicted for e = 0.01 in Figure 
[S] The corresponding Poincare section is shown in 
Figure (The 4 : 1 resonance bands are not shown 
in this plot.) 

The predictions for centers are not as good as those 
for saddles, but there is nevertheless good quantita- 
tive agreement between observation and prediction, 
especially considering that a leading-order perturba- 
tion method has been employed. Of course, given 
that higher-order resonances occupy smaller regions 
of phase space, the absolute errors indicate that these 
predictions are not as good as the same absolute er- 
rors would be when studying lower-order resonances. 
This caveat notwithstanding, our theoretical analysis 
does an excellent job of determining the location of 
resonances and offers a useful tool for locating high- 
order resonances (and thus studying band structure 



in great detail) in numerical simulations. 

7 Conclusions 

In this paper, we studied in depth the band struc- 
ture of BECs in periodic lattices. We approached 
this problem using a coherent structure ansatz, in 
contrast to the Bloch wave ansatz of earlier stud- 
ies. 26 - 39 ' 52 

Using a technically delicate perturbative approach 
relying on elliptic function solutions of the integrable 
NLS, we examined the spatial resonance structure 
(band structure) of coherent structure solutions of 
the NLS in considerable detail, providing both an 
analytical description and numerical verifications of 
this theory. We derived conditions for the onset of 
2m' : 1 spatial resonances for all integer m' and de- 
veloped analytical expressions for the width of these 
resonance bands and the locations of saddles and cen- 
ters therein. Comparison with numerical simulations 
of primary, secondary, and tertiary resonances illus- 
trate the applicability of our analytical theory. 

Utilizing a simpler perturbative approach that em- 
ploys Lindstedt's method and multiple scale analysis, 
we also established wave number-amplitude relations 
for coherent structure solutions of the NLS with a pe- 
riodic potential. In so doing, we explored 2 : 1 spatial 
resonances and illustrated the utility of phase space 
analysis for the study of band structure as well as the 
structure of modulated amplitude waves in BECs. 

In sum, our perturbative approach does an excel- 
lent job of determining the location of resonances 
and analyzing their structure and offers a useful tool 
for locating high-order resonances (and thus study- 
ing BEC band structure in great detail) in numeri- 
cal simulations. An important open direction, to be 
addressed in a future publication, is the extent to 
which the theory developed here is an effective start- 
ing point for studies of the chaotic dynamics of BECs. 
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Equilibrium at (0, 0) 


Equilibria at (i2*,0) ^ (0,0) 


9 


fj, 


Center 


None 




+ 


Center 


None 





+ 


Center 


Saddles 


+ 


+ 


Saddle 


None 


+ 




Saddle 


None 







Saddle 


Centers 







Table 1: Type of equilibria of (JSJ) when V(x) — 0, and c = 0. 



e = 0.01 e = 0.05 



Quantity 


Perturbative 


Numerical 


Perturbative 


Numerical 


Y x 


0.28133 


• 


0.28133 




Y s 


0.27949 


• 


0.27213 




R s 


±0.7477 


±0.66 


±0.73775 


±(0.56- 0.66) 


Y c 


0.28317 


• 


0.29053 




R c 


±0.7526 


±0.757 


±0.76227 


±0.774 


Y 


0.16175 


• 






Y 

1 max 


0.40091 


• 






Rin 


±0.56877 


±0.66 






Rout 


±0.89545 


±0.85 






Ymin,2 


0.22203 


• 






Y max ,2 


0.35807 


• 






Rin, 2 


±0.66638 


±0.66 






Rout,2 


±0.84626 


±0.85 







Table 2: Comparison of perturbation theory and numerics for 2 : 1 resonances (k = 1.5, i5 = 1). In this 
Table, Y\ is the action value of the primary resonance, Y s is the location of its nearby saddle, (R s ,0) is its 
location in (R, S)-coordinates, Y c is the location of the nearby center, (R c , 0) is its location in (R, 5')-space, 
Y min is the minimum action value of the separatrix determined from l|81f) . Y max is the maximum determined 
from Q82f) . Ri n is where the inner separatrix crosses the i?-axis, R ou t is where the outer separatrix crosses 
the i?-axis, Y m i n ^ and Y max .2 are the minimum and maximum actions obtained by solving (|77|l numerically, 
and Rin,2 and R ou t.2 are their corresponding predictions of where the inner and outer separatrices cross the 
i?-axis. The symbol • means a calculation is not applicable and — means it was not computed. 
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e = 0.01 e = 0.05 



Quantity 


Perturbative 


Numerical 


Perturbative 


Numerical 


Yi 


0.19443 


• 


0.19443 


• 


Y s 


0.18499 


• 


0.14718 




R s 


±27202 


±(0.19-0.20) 


±0.24264 




Y c 


0.20389 


• 


0.24169 




R c 


±0.28558 


±0.2908 


±0.31093 


±0.335 


Y 

1 min 


0.07485 


• 






Y 

1 max 


0.31402 


• 






Rin 


±0.17303 


±0.19 






Rout 


±0.35441 


±0.37 






Y m i n ^2 


0.09644 


• 






Y max ,2 


0.34904 


• 






Rin, 2 


±0.19640 


±0.19 






Rout,2 


±0.37366 


±0.37 







Table 3: Comparison of perturbation theory and numerics for 2 : 1 resonances (k = 0.75, 5 = 0.2). The 
quantities computed are defined in the caption of Table [21 



e = 0.01 e = 0.05 



Quantity 


Perturbative 


Numerical 


Perturbative 


Numerical 


Y 2 


0.37358 


• 


0.37358 


• 


Y s 


0.37294 


• 


0.37036 


• 


R s 


±0.86364 


±0.88 


±0.86065 


±0.88 


Ss 


±0.68389 


±0.687 


±0.68293 


±0.68 


Y c 


0.37422 


• 


0.37680 


• 


(R C ,S C ) 


See text 


(±0.691, ±0.332) 


See text 


(±0.697, ±0.330) 


Y 

1 min 


0.34814 


• 


0.31670 


• 


Y 

1 max 


0.39902 


• 


0.43046 


• 


Ymin,2 


0.35571 


• 


0.32989 


• 


Y max ,2 


0.41237 


• 


0.46240 


• 



Table 4: Comparison of perturbation theory and numerics for 4: 1 resonances (k = 2.5, 6 = 1). Tha action 
value of the secondary resonance is denoted Y2 • Saddles that intersect the i?-axis are denoted (R s , 0) , and 
those that intersect the S'-axis are denoted (0,S S ). Centers are denoted (R C ,S C ). The other quantities 
computed are defined in the caption of Table [5] 



Quantity 


Perturbative 


Numerical 


Y 3 


0.36857 


• 


Y s 


0.36851 


• 


Rs 


±0.85850 


±(0.859- 0.860) 


Y c 


0.36863 


• 


Rc 


±0.85864 


±0.870 


Y 

1 min 


0.36214 


• 


Y 

1 max 


0.37500 


• 


Ymin,2 


0.36614 


• 


Yfnax.2 


0.38653 


• 



Table 5: Comparison of perturbation theory and numerics for 6 : 1 resonances (k = 3.8, 5 = 1, e = 0.01). 
Tha action value of the tertiary resonance is denoted I3. The other quantities computed are defined in the 
caption of Table [21 
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-1.5 1.5 
R 

(c) 

Figure 1: Phase portraits of coherent structures in BECs with no external potential. The signs of ^ and g 
determine the dynamics of (0. (a) Repulsive BEC with /i > 0. The two-body scattering length is a = 0.072 
nm, the value 28 for atomic hydrogen ( 1 H). Orbits inside the separatrix (which consists of two heteroclinic 
orbits) have bounded amplitude R(x). The period of such orbits increases as one approaches the separatrix, 
whose period is infinite, (b) Attractive BEC with [i > 0. The two-body scattering length is a = —0.9 nm, 
the value 25 ' 47 for 85 Kb. (c) Attractive BEC (again 85 Kb) with /i < 0. Here there are two separatrices, each 
of which encloses periodic orbits satisfying R ^ 0. 
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Figure 2: As the wave number k of the perturbation is increased, periodic behavior persists for larger 
\Vq\. The initial condition in this plot is (i?(0), 5(0)) = (0.05,0.05), and the parameter values a — 0.072, 
fi = 10, m = 0.5, and xq = are used for each trajectory, (a) Poincare section determined by sin(ra) = 0. 
Trajectory (1) corresponds to (k, Vq) — (100, 10), trajectory (2) to (k, Vq) = (100, 100), and trajectory (3) to 
(k, Vq) = (10, 10). These quasiperiodic solutions indicate the existence of nearby periodic orbits, (b) Phase 
space plots of the trajectories in (a). Trajectory (1) is the closest to being periodic and trajectory (3) is 
the furthest away, (c) Amplitude R as a function of space x for trajectories (l)-(3). The band structure 
of BECs can be studied not only in real space but also in phase space by plotting Poincare sections and 
trajectories, as indicated in (a) and (b). Examining the proximity of a trajectory to periodicity is most 
easily accomplished in phase space, (d) Coherent structure corresponding to quasiperiodic trajectory (1). 
This plot depicts Re(ip). The horizontal axis represents time, and the vertical one represents space. The 
darkest portions are the most negative, and the lightest are the most positive. 



20 




Figure 3: Poincare section for the parameter values fi = —10, m = 0.5, xq = 0, Vq = 5, k — 10, and 
a = —0.9 nm, corresponding to the experimentally determined scattering length 25 ' 47 for 85 Rb. The depicted 
trajectories include examples which are quasiperiodic, locally chaotic (near the resonances), and globally 
chaotic (the stochastic sea). 




Figure 4: (a) Spatial profile of the coherent structure corresponding to the locally chaotic trajectory in 
FigureEl The initial conditions are (i?(0), 5(0)) w (-0.01818215,-5.23268358). (b) Spatial profile of the 
coherent structure corresponding to the globally chaotic trajectory in Figure |3J The initial conditions are 
(i?(0),5(0)) w (-1.13283530,1.28334013). 
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(c) _3 _2 -1 12 3 (d) -3 -2 -1 12 3 



Figure 5: Poincare sections (a), (b) and resonance Hamiltonians K\ (c), (d) for k — 1.5 and 5 = 1. (a) 
Poincare section for e = 0.01. The 2 : 1 resonances are displayed, as indicated by the numbered trajectories, 
(b) e = 0.05. (c) Resonance Hamiltonian for e = 0.01 with vertical axis in units of action Y and horizontal 
axis in units of £/J'(Y). (d) e = 0.05. 
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-3-2-10 1 2 3 
Figure 6: Resonance Hamiltonian K2 for k = 2.5, 6=1, and e = 0.05. 




-3-2-10 1 2 3 
Figure 8: Resonance Hamiltonian K3 for k = 3.8, 6 = 1, and e = 0.01. 




0.22 



(a) 



-0.22 




0.81 
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Figure 9: (a) Poincare section for k — 3.8, 5=1, and e = 0.01. (b) Close-up of the resonances in (a). Both 
6:1 (1) and 8:1 (2) resonances are displayed. A higher-order resonance (3) is also depicted. Although not 
shown, 4 : 1 resonances are also present for this choice of parameter values. 
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